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The ground state of Bose-Einstein condensates with attractive particle interaction is metastable. 
One of the decay mechanisms of the condensate is a collapse by macroscopic quantum tunneling, 
which can be described by the bounce trajectory as solution of the time-dependent Gross-Pitaevskii 
equation in imaginary time. For condensates with an electromagnetically induced gravity-like in- 
teraction the bounce trajectory is computed with an extended variational approach using coupled 
Gaussian functions and simulated numerically exact within the mean-field approach on a space-time 
lattice. It is shown that the variational computations converge very rapidly to the numerically exact 
result with increasing number of Gaussians. The tunneling rate of the condensate is obtained from 
the classical action and additional parameters of the bounce trajectory. The converged variational 
and numerically exact results drastically improve by several orders of magnitude the decay rates 
obtained previously with a simple variational approach using a single Gaussian-type orbital for the 
condensate wave function. 

PACS numbers: 03.75.Hh, 34.20. Cf, 34.80. Qb, 04.40. b 



I. INTRODUCTION 

The macroscopic occupation of the bosonic ground 
state of ultracold quantum gases has been predicted by 
Bose and Einstein. At least since the first realization 
of Bose-Einstein condensates (BECs) in 1995 they 
are a central part of experimental and theoretical atomic 
physics. In a cold diluted Bose gas the particles interact 
via short-range contact interactions determined by the 
s-wave scattering length a, which can be experimentally 
varied with the help of Feshbach resonances. In addition 
long-range interactions exist, e.g., in dipolar condensates 
[3], or in condensates with an attractive l/r interaction, 
which have been proposed to be electromagnetically in- 
duced with a setup of six laser triads ^ . The condensates 
are typically held in a harmonic trap, however, in case of 
an attractive l/r interaction the BEG can be self-trapped 
[316], which means that no external trap is needed. 

When the scattering length of the contact interaction 
is varied, the stable ground state plus an unstable ex- 
cited state of a BEG with attractive l/r interaction are 
created in a tangent bifurcation at a critical scattering 
length a = Ocr [3 E] ■ Below the critical scattering length 
the contact interaction is so attractive that the conden- 
sate collapses and no stationary state does exist [8 . The 
collapse of the condensate, e.g., induced by thermal exci- 
tations [9HT2] , is also possible at scattering lengths above 
the bifurcation point. However, even without thermal 
excitations, i.e., at zero temperature T = 0, the ground 
state of the condensate is only a metastable state and 
can collapse by macroscopic quantum tunneling as long 
as the contact interaction is attractive. The process can 
be described by using the bounce trajectory as outlined 
by Stoof [13] and by Freire and Arovas [T4] . 

The aim of the present paper is a thorough investi- 
gation and computation of the bounce trajectory of a 
BEG. It has been shown that the stationary states of 



the GPE can be excellently described with an extended 
variational approach using coupled Gaussian functions 
[T51 [TC] . Here we will demonstrate that computations 
with coupled Gaussians are also a full-fledged alterna- 
tive to simulations on a space-time lattice to obtain the 
bounce trajectory which describes the macroscopic quan- 
tum tunneling of the BEG. We apply the formalism to a 
self-trapped BEG with an attractive l/r interaction. On 
the one hand the symplicity of this spherically symmetric 
system allows for an intuitive access to the bounce tra- 
jectory similar to the one-dimensional potential picture 
used by Stoof [l^ and Freire et al. [MKIT]. On the other 
hand we can show that it is possible to include long-range 
interactions into the instanton formalism for BEGs. The 
theory presented in this article can be applied to more 
realistic dipolar condensates 4J as well since their de- 
scription does not contain a conceptional difference. 

The paper is organized as follows. The basic concepts 
are briefly outlined in Sec. [H] In Sec. |III A the time- 
dependent variational principle (TDVP) is applied to an 
ansatz of coupled Gaussians for the wave functions in 
imaginary time. Periodic solutions of the equations of 
motion for the variational parameters are obtained by 
application of a multi-shooting algorithm, which con- 
verge with increasing period to the bounce trajectory. 
In Sec, nil B] the method for exact numerical simulations 
on space-time lattices is elaborated. The results of both 
methods are presented and compared in Sec. jlVj G on- 
cluding remarks and an outlook are given in Sec. |V| 



II. BASIC CONCEPTS 

For condensates without long-range interactions the 
macroscopic quantum tunneling has been investigated 
by Stoof jl3j using a variational approach with a single 
Gaussian function as a simple approximation. For a self- 
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trapped BEC with an attractive l/r interaction instead 
of an external harmonic potential the ideas outlined in 
Ref. [13] can briefly be summarized as follows. Starting 
from the many-body Schrodinger equation a mean-field 
approach leads to the time-dependent nonlinear Gross- 
Pitaevskii equation (GPE) 



i-V(r, t) = [-A + V,{r, t) + K,(r, t)] ^{r, t) (1) 
with the contact and long-range potentials 



K(r,t) = 87r7Va|V(r,t)r , 
K(r,t) = -2iV /dV''^^"''')'' 



r — r' 



(2) 
(3) 



where a is the scattering length, N is the number of 
bosons, and "natural" units 6J were used. Lengths are 
measured in units of a„ = h? /{mu), energies in units of 
Eu — h? /(2ma^), and times in units oit^ = h/Eu, where 
u determines the strength of the atom-atom interaction 
[5] , and m is the mass of one boson. Exploiting the scal- 
ing properties presented in Ref. ^ and introducing the 
scaled variables 



(r, a, f, V-) = {Nr, N'^a, NH, N~^/^iP) 



(4) 



the GPE keeps the form given in Eq. ([T]), however, with 
the particle number in the potentials Vc and in Eqs. ^ 
and ^ formally set to iV = 1. Note that the scaled GPE 
for the self-trapped BEC only depends on one parameter, 
viz. the scaled scattering length. In the following we use 
the scaled variables unless otherwise stated and omit the 
tilde. 

Using a Gaussian-type orbital 
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normalized to jj^/'lp = 1 as a simple variational ansatz for 
the condensate wave function the TDVP can be applied, 
and it can be shown [5j that the parameters q and p obey 
the canonical equations of motion 



dH 
dp 



dH 



(6) 



with the Hamiltonian H{q,p) = + V{q) and the po- 
tential 
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Important features of the BEC can be simply obtained, 
at least qualitatively, from the potential V{q) shown in 
Fig. [IJa). For scattering lengths a < acr = — 37r/8 — 
— 1.1781 the potential V{q) does not possess any station- 
ary points, i.e., the condensate collapses. For < a < 
two stationary points exist, the local minimum at q — 
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tunnel process , 
, '^9''V(q) 



FIG. 1. (Color online) (a) Potential V{q) of a self-trapped 
BEC with attractive l/r interaction at various values of the 
scattering length a. (b) Scheme of the inversion of the poten- 
tial and the periodic orbits approaching the bounce trajectory. 



Qmin and maximum at q = (/max of the potential indi- 
cate the stable ground state and the unstable excited 
state, with V{qinin) and l^((?niax) the mean-field energies 



Emf - 

lively. 



h{Vc + Vu)\'tp) of the two states. 



respec- 



From Fig.[TJa) it is obvious that the potential V{q) at 
small values of q drops below the mean-field energy of 
the ground state. The path through the potential bar- 
rier connecting the position of the ground state and the 
collapsing region [see the arrow in Fig. [ija)] is classically 
forbidden. In quantum mechanics, however, the barrier 
can be crossed by a tunneling process, which is described 
by propagation along that path, called the bounce trajec- 
tory, in imaginary time. The decay rate of the metastable 
ground state depends on parameters of the bounce tra- 
jectory. 

In the mean-field approach the many-body ground 
state is the product of N identical one-particle wave func- 
tions, and thus the collapse of a BEC requires the col- 
lective tunneling of all N bosons. As a consequence for 
macroscopic quantum tunneling an effective Planck con- 
stant h/N must be used in the formula for the decay rate 
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Tq. The final result derived by Stoof [TH] reads 



where 



5b = 2 / ^Viq) - Viqnun)dq 



(8) 



(9) 



is the Euclidean action of the bounce trajectory (with 
the turning points gb and Qmin) in imaginary time it — > t 
running in the inverted potential i.e., V{q) — >■ —V{q) 
as illustrated in Fig. [ijb). The parameter cjq is the fre- 
quency of oscillations in the potential minimum around 
q ~ gmin given by 



Wo = — V^"(gmin) 



(10) 



and Wo is defined by the condition that the bounce trajec- 
tory approaches the maximum of the inverted potential 
—V{q) [see Fig. [l]^b)] in the hmit r — >■ ±oo as 



ujo 



(11) 



The computation of the decay rate Fq in Eq. Q with 
the parameters S'b, ujq, and vq of the bounce trajectory 
as explained above is only a rough approximation for 
the following two reasons. Firstly, the ansatz with a 
single Gaussian function in Eq. ([s]) cannot describe the 
true condensate wave function. The variational solutions 
significantly deviate from those found in full numerical 
simulations even for the stationary ground and excited 
state 0[7]. Secondly, the true solution of the GPE Q 
can fluctuate around both the stationary ground state 
and the time-dependent bounce trajectory. The modes 
of these fluctuations are described, in leading order, by 
the Bogoliubov-de Gennes equations [18] . In Eq. ([s]) only 
fluctuations along the direction of the path are consid- 
ered, i.e., for a system with only one degree of freedom. 
The true condensate wave function has many degrees of 
freedom, which allows for fluctuations perpendicular to 
the direction of the path. The fluctuation determinant 
resulting from those additional degrees of freedom is not 
accounted for by Eq. ([s]). 

The method proposed by Stoof [13] requires the de- 
scription of the bounce trajectory with canonical vari- 
ables and the Hamiltonian formalism. Such a descrip- 
tion is, however, restricted to a variational ansatz for the 
condensate wave function with a single Gaussian function 
as in Eq. ([5]), and thus any improvement of the simple 
variational approach is a nontrivial task. 

A decisive step for the improvement of the bounce tra- 
jectory approach was achieved by Freire and Arovas |14j . 
They found a way to compute the bounce process beyond 
the interpretation in terms of a classical trajectory by ap- 
plying the instanton formalism to a BEG in a harmonic 
trap without long-range interaction. The bounce trajec- 
tory was simulated numerically on a space-time lattice 
without any restriction to the wave function. 



To solve the GPE ([!]) in imaginary time it — )■ r the two 
fields tp{r, t) and ip{r, r) are introduced, which obey the 
equations 



dr 
with 



H{r,T) 
K(r,r) 



H^Pir,T), —i}{r,T) = H^{r,T), (12) 



-A + V;(r,T) + K(r,r)-M, (13) 
87riVaV;(r,T)-0(r,T) , (14) 



-2N 



r — r' 



(15) 



Note that /i is a free parameter in Eq. ( 13 ) and the func- 



tion %p{r^T) in imaginary time plays the role of the com- 
plex conjugate tp* (r, t) in real time dynamics. The prob- 
ability density is given by ijiil) and the norm defined by 
/ V'(r, t)'0('", ''■)d^r is conserved. The bounce trajectory 
(with real valued ■0 and ■)/') must fulfill the boundary con- 
ditions [Hdl] 



(16) 



lim ^{r,T)= lim tp{r,T) = il;g^{r), (17) 

r— f±oo r— f±oo 

with '4'g{r) the wave function of the stationary ground 
state. The Euclidean action of a periodic trajectory in 
imaginary time with period /? is given by 



13/2 



fS/2 



dr / dV { V'— - n 



with the associated Hamiltonian density 



H = -VV'VV' + 47ra (V'V')^ ~ / ^ 



3 ,1p'lp'tj)ll) 

r 

\r — r'\ 



(18) 



(19) 



Due to the symmetry ( 16 1 oitp and if) and the behavior of 
5c [qh] and [qo] the difference A5o = [q^] — [go] be- 
tween the actions of the bounce trajectory and the fixed 
point is obtained as [14 



A5c = 



13/2 



dr / (5?r < xp—^ — ih—^ 
dr dr 



(20) 



Within the mean-field approach the formalism allows, in 
principle, for the exact computation of the bounce tra- 
jectory and related parameters needed as input for the 
calculation of the decay rate of the BEG in Eq. Q with- 
out resort to the ansatz ([5| and any canonical variables. 



III. CALCULATION OF THE BOUNCE 
TRAJECTORY 

The computation of the exact bounce trajectory which 
describes within the mean-field approach the macroscopic 
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quantum tunneling of a metastable BEC is the central is- 



sue of this paper. We first introduce in Sec. Ill A the ex- 
tended variational approach with coupled Gaussian func- 
tions, and then in Sec. |IIIB] the numerical simulations on 
space-time lattices. The derivations are presented for a 
self-trapped BEC with attractive 1/r interaction, where 
the condensate wave function is spherically symmetric. 



A. Variational approach with coupled Gaussians 

The variational ansatz with a single Gaussian func- 
tion in Eq. ([5| cannot exactly reproduce the true conden- 
sate wave function. However, it has been shown that an 
extended variational approach with coupled Gaussians 
yields excellent results for both the stationary ground 
and excited state of the BEC [TB] . Only very few (about 
three to five) coupled Gaussians are sufficient to achieve 
convergence of the wave function, the mean-field energy, 
and the chemical potential to those obtained in numeri- 
cal computations on grids. We therefore will use a varia- 
tional approach with coupled Gaussians also for the cal- 
culation of the bounce trajectory. For the dynamical cal- 
culations in imaginary time we write both functions ip 
and ip ^ superpositions of Ng Gaussian functions, 



we obtain the condition 

51^ [ d^r {S^ [(f> + HiP] +S(f>[4>- H^p] } 

d'ip 



/d3.^ 

J, 1 



fc=l 



(0 -I- H^) 

9-0 



Szk 



dzk 



dzk_ 

Szkl^O, (25) 



and thus, because all Sz^ and Szk are independent 

dip 



dV 



j d^r {ij - HiP^ 



dip 
dzk 



(26a) 
(26b) 



for fc = 1, . . . ,2Ng. Using the ansatz (21) with coupled 



Gaussians and the Hamiltonian (13) in Eq. (26 1 a calcu 



lation similar to that in Ref. |15j vields the equations of 
motion for the variational parameters 



(2) 



7fc = 6Afc + - ^J. 



Ak = ^Ai - 
7fc = -&Ak 



-(2) 



(27a) 
(27b) 
(27c) 
(27d) 



i'ir, r) = ^ e +')''= M] = ^ gfe , (21a) ^i^j^ fc = 1, . . . , TV^. As in Eq. (p| ^ is a free parame- 



fc=l 



ter, which becomes the chemical potential for stationary 



^(r, t) = ^ e~[Ak(-r)r''+lk(T)] = ^ . (21b) 

fc=l 



states. The parameters u^^'', v^^\ v'^^^ and v'"^' in Eq. ( |27[ ) 
are obtained by solving the two linear systems of equa- 
tions 



,(0) 



fc=i 



The ANg time-dependent variational parameters can be 
written as two vectors 



(o)rii , Y^„,(2)r^2i 



z 
z 



(^i,...,^w^,7i,...,7wj , (22a) 
(ii,...,i« ,7i,...,7« ) . (22b) 



fc=l k=l 



Ik — 



Ik 



fc=i 



k=l 



k=l 



The equations of motion for the variational parameters z 
and z can be obtained using the TDVP of Dirac, Frenkel, 
and McLachlan [19l [20], which must be modified and 
adopted in a special way for the dynamics in imaginary 
time. In imaginary time the functional 



and 



N, N, 

n(")rii , Y^„-i(2)rJ 



fc=i 



N„ 



fc=l 

Na 



k=l 
No 



I = 



dV [<P - HiP] [(j) + HiP] 



(23) 



must be minimized with respect to (p = ip and (p = ip, 



with the Hamiltonian H given in Eq. ( 13 ). The dots now 
mark derivatives d/dr. With 



E 

k=l 



dip 
dzh 



Szi 



2Ng 

E 

k=l 



dip 
dzk 



6zi 



(24) 



(28a) 
(28b) 

(29a) 
(29b) 

(30) 



All required integrals can be calculated analytically and 
are given in Appendix [A] 

The equations of motion ([27]) are a set of ANg or- 
dinary differential equations, which can be integrated. 



E4"^w- + E<'^[-']- = E[^]-' 

fc=i 

N, 



k=l k=l 



fc=l 



with V ^ Vc + Vu and the notation 

[0]ik = {9i\0\gk) = I d\-giOgk 
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e.g., with a Runge-Kutta method. Because the period 
/3 of the bounce trajectory is infinite we start search- 
ing for periodic orbits, with the boundary conditions 
'^{t = 0) = i^ir = 0) and ^{t = /3/2) = ^{t = /3/2), 
near the excited state, as illustrated in Fig. [ijb), and 
then increase /3. In the limit /3 — > oo the periodic or- 
bits converge towards the bounce trajectory. The pe- 
riodic orbit search could be performed, in principle, by 
varying the initial conditions of the variational param- 
eters at r = with ip{0) — V'(O) in such a way, that 
■0(/3/2) — ■!/)(/3/2) is fulfilled at the other turning point. 
Numerically, however, this is an impossible task because 
the equations of motion are extremely unstable. The rea- 
son is that, e.g., for the ground state all stable modes of 
the Bogoliubov spectrum in real time become unstable 
modes in imaginary time, and thus all perturbations of 
the ground state (or a trajectory) increase exponentially 
in imaginary time. To overcome this problem we use a 
multi-shooting algorithm, where the trajectory is divided 
into small segments. Using M segments the total set of 
AMNg + 1 parameters [including the parameter /i in the 



equations of motion (27)] can be determined to fulfill the 



AMNg continuity conditions of the variational parame- 
ters and the condition J d^ripip = 1 for normalization 
by using a multidimensional Newton's method. Because 
the norm of 'tp^p is conserved, the total set of parameters 
can be slightly reduced by introducing the parameters 
Ik = 7/c - 71 and % = 7^ - 71 for fc = 2, . . . , Ng instead 
of 7fe and 7^ for A; = 1, . . . , Ng. 

As an example Fig. [2] shows the variational parame- 
ters of a periodic trajectory at scaled scattering length 
a = —0.9 with period /3 = 47.4 near the bounce. The 
BEC is described by a superposition of A^g = 3 Gaussian 
functions. The high values of the parameters Aj^ and A^. 
in Fig.[2|a) indicate a small width of the condensate wave 
function at the bounce around t = 0. 

Once a bounce trajectory has been found improved val- 
ues of the parameters S'b, wqi and Vq in the formula ([s]) 
for the decay rate of the BEC can be determined. The 
Euclidean action S'b = A.S0 of a periodic trajectory is 
given in Eq. (20 1. For the ansatz with coupled Gaus- 



sians the action can directly be expressed in terms of the 
variational parameters, 

AS, = dTj2 

X [(^Ai-Ak){gi\r^\gk) + iii-ik){gi\gk)} ■ (31) 



The matrix elements are given in Eqs. (All and (A2| in 
Appendix [A] 

The frequency ujq is obtained, instead of using the po- 
tential in Eq. (10 1, by a stability analysis of the ground 
state. For the ansatz with coupled Gaussian functions 
the stability parameters are obtained as eigenvalues of 
the Jacobian matrix constructed by varying the varia- 
tional parameters around the fixed point solution |15j . 
The eigenvalue ujq belongs to the variation of the param- 
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FIG. 2. (Color online) Variational parameters (a) ^^(r), 
^fc(r) for fc = 1, . . . , Ng and (b) 7fc('r)-7i('r). ^k{T)-ji{T) for 
fc = 2, . . . , Ng of a periodic trajectory with period /3 — 47.4 
at scaled scattering length a = —0.9 described by Ng — 3 
coupled Gaussian functions. 



eters in the direction of the bounce trajectory. 



The parameter vq cannot be determined using Eq. ( 11 ) 



because no canonical coordinate g(r) exists for the ansatz 
with coupled Gaussians. However, when the bounce tra- 
jectory is approached via periodic orbits, vq can be esti- 
mated using the mean- field energies and E^f{P) of 
the ground state and periodic trajectories with period /3. 
Using a model with an inverted harmonic potential for 
approaching the saddle in imaginary time we obtain (see 
Appendix Ib|) 



Vq — lim 



BMP) - 
2m 



."0/3/2 



(32) 



Results obtained with coupled Gaussians are presented in 
Sec. |IV| and compared with numerically exact simulations 
on space-time lattices. 



B. Numerical simulations on space-time lattices 

To verify the validity of the extended variational ap- 
proach with coupled Gaussian functions we have carried 
out numerically exact computations by simulating the 
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bounce trajectory on a space-time lattice. Because the 
condensate wave function is spherically symmetric only 
the r coordinate needs to be discretized in space. We 
set ipij = tp{iAr,jAT) and tpij = iJj{iAr, jAr) with 
i = 1, . . . ,m and j = 1, . . . , n for an (m, n) space-time 
lattice with step sizes Ar and At, respectively. 

The propagation in imaginary time is obtained with 
the operators 



U 



-H7 



U 



Hi 



(33) 



applied to the wave functions ip and V', respectively. The 
Hamiltonian 7J = T + Vc -f Vu — /i is given in Eq. ( 13 1. 
Because the operators of the kinetic energy T = —A and 
of the potential V = Vc-\-V^ do not commute we adopt 
the Baker-Hausdorff formula 



-{T+V)At 



„-TAt -V\t/2 



(34) 



valid for small time steps At and apply the split-operator 
method. The operators q^^^'^/'^ are applied in coordi- 
nate space, and the operators e^"^'^'^ in momentum space. 
The Fourier transform of a spherically symmetric wave 
function tp{r) yields 



An 

T 



riplr) sin(A:r)dr , (35) 



which can be computed efficiently on the grid by a fast 
sine transform |2lj. For the potential Vu(r) of the attrac- 
tive 1/r interaction [see Eq. (15l] we obtain by applica- 



tion of the convolution theorem 



K(0 = -^"M HH^J" 



1 

fc2 



r'4>{r')il^{r') sin(fcr')dr' 



sin(fcr)dfc , (36) 



where the integrals are also computed by fast sine trans- 
forms. For large values of r the potential Vu{r) asymp- 
totically approaches — 2/r, however, the sine transform 
yields Vu(?'max) = at the border r^ax ~ mAr of the 
grid. The error is corrected by a shift of the potential 
K(r-)^K(r)-2/w. 

The space-time lattice can be initialized, e.g., by dis- 
cretizing a trajectory obtained with the extended varia- 
tional approach introduced in Sec. |III A| The remaining 
task is then to fulfill the matching conditions that the 
grid points of ip and ip at time t = jAt propagated by 
a time step At coincide with the grid points of ip and 
■0 at time t = (j + 1)At. For a periodic trajectory fur- 
thermore the boundary conditions at the turning points 
'(/'(t = 0) = iplr ~ 0) and = ?^At) = ^(t — nAr) 
must be fulfilled. The matching conditions yield a high- 
dimensional system of 2m{n — 1) coupled equations, 
which can be solved by a Newton method. For a typical 
size of the space-time lattice with to = 64 and n = 151 
used in our calculations the Jacobian matrix in Newton's 
method has the dimension 19200 x 19200, which means 
about 3.7 X 10* matrix elements. However, the matching 




FIG. 3. (Color online) Density [ijjip]{r) of a condensate wave 
function for a periodic trajectory with period (5 = 42.45 com- 
puted on the space-time lattice at increasing values |t| (from 
top to bottom line) of the imaginary time. The scaled scat- 
tering length is a = —0.9. 



conditions only depend on the grid points at two adjacent 
time steps, and thus the Jacobian has a well pronounced 
band structure which allows for an efficient solution of the 
linear system of equations |21j . Furthermore, the norm 
conservation can be used to slightly reduce the number 
of parameters. 

Similar as discussed in Sec. IIII Al for the variational 
calculations we start searching for periodic trajectories 
on the space-time lattice with a short period /3 at mean- 
field energies slightly below the excited state (see Fig. [T]) 
and then increase /3. As an example Fig. [3] shows the 
density ['iptp]{r) of a condensate wave function for a pe- 
riodic trajectory with period /3 — 42.45 computed on 
the space-time lattice at various values of the imaginary 
time T. The wave function is narrowest at the bounce at 
T = and becomes broader with increasing time. Note 
that the area under each curve in Fig.|3]is not the square 
norm 27r r'^ip(r)tp{r)dr and thus not conserved. As 
discussed above, the bounce trajectory is obtained in the 
limit /3 — >■ oo. 

The Euclidean action of a periodic trajectory on the 
space-time lattice is computed with Eq. (20 1. The ex- 



act value of the frequency ujq is obtained by numerically 
solving the Bogoliubov-de Gennes equations at the sta- 
tionary ground state as discussed in Refs. [8l[22]. The 
parameter vq is determined with Eq. ( 32 1 . 



IV. RESULTS AND DISCUSSION 

We now present and compare the results obtained 
with the extended variational approach using coupled 
Gaussian functions and by the numerically exact sim- 
ulations on space-time lattices. As outlined above the 
bounce trajectory is approached via periodic trajecto- 
ries fulfilling the boundary conditions ipifi) = '0(0) and 
iIj{/3/2) = ■i/;(/3/2) at the turning points in the limit of 
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FIG. 4. (Color online) (a) Euclidean action of periodic tra- 
jectories at constant scaled scattering length a = —0.9. (b) 
Euclidean action of the bounce trajectory as function of the 
scaled scattering length a. Shown are the results obtained by 
the variational approach using up to = 5 coupled Gaus- 
sians and by numerically exact simulations on grids. The 
variational results with Ng > 3 agree within the line width 
with the grid computations. 



FIG. 5. (Color online) Decay rate by macroscopic quantum 
tunneling of a BEC with (a) iV = 30 and (b) iV = 100 par- 
ticles obtained by Eq. (|8| with the parameters S'b, ojo, and 
uo of the bounce trajectory. The variational calculations with 
Ng = 1 to 5 coupled Gaussians converge rapidly to the grid 
computations, for Ng > 3 deviations are less than the line 
width. 



an infinite period /3 — >■ oo. The periodic orbit search is 
started close to the stationary excited state as illustrated 
in Fig. [l|b) where the period /3 = 2tt/u)c is given by the 
corresponding stability eigenvalue Wc of that state. The 
periodic orbit is then varied towards the bounce trajec- 
tory by increasing j5 in small steps. 

The Euclidean action of the periodic orbits at constant 
scaled scattering length a ~ —0.9 as function of the pe- 
riod 13 is shown in Fig. |4]ja). The results are obtained 
by the variational approach using up to A'g = 5 cou- 
pled Gaussians and by numerically exact simulations on 
space-time lattices. Evidently, the variational calcula- 
tions converge rapidly to the solution of the numerical 
grid calculations for an increasing number of Gaussians. 
For more than three Gaussians the deviations are less 
than the line widths. However, it is important to note 
that the converged action deviates significantly from the 
result obtained with a single Gaussian function [that ac- 
tion is ^ 0.83 for a = —0.9 and thus far outside the 
region presented in Fig. |4][a)] . The action decreases with 
increasing period, but becomes nearly constant in the re- 



gion /3 > 20, which indicates the rapid convergence to 
the action of the bounce trajectory. 

The limit /3 — oo of the Euclidean action has been 
estimated for periodic trajectories at various scaled scat- 
tering lengths. The resulting action of the bounce tra- 
jectory is presented in Fig. [4jb). Again, the variational 
results converge rapidly to those obtained by simulations 
on space-time lattices. The action of the bounce trajec- 
tory vanishes at the critical scattering length, where the 
ground and excited state emerge in a tangent bifurcation. 
The exact critical scattering length is = —1.0251, 
however, in the variational calculation with a single 
Gaussian the bifurcation point is significantly shifted 
to Act'"^ — —1.1781. As a consequence in the region 
a > a°f = —1.0251 the difference between the action ob- 
tained with Ng ~ 1 and the exact action is ^ 0.5, i.e., 
quite large. 

The Euclidean action shown in Fig. |4|^b) and the pa- 
rameters Wo and vq of the bounce trajectory can now be 
inserted into Eq. ([s]) to determine the decay rate Fq of the 
BEC with attractive l/r interaction. Although particle 
number scaled units [see Eq. ^] are used in Eq. ([8| the 
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TABLE I. Decay rates To in SI units for self-trapped conden- 
sates of A'^ = 30 and A'' — 100 *^Rb atoms with the parame- 
ters [6l|23| a„ = h^/{mu) ~ 2.5 x 10"" m and tn = 2mal/h ~ 
27.1s. 





p{Ar.30)j^-lj 




1 


3.5 X 10"'^ 


9.65 X 10"^" 


2 


0.70 


0.066 


3 


1.18 


0.45 


4 


1.21 


0.49 


5 


1.21 


0.50 


grid 


1.21 


0.50 



decay rate depends, through the replacement h ^ h/N 
for the macroscopic quantum tunneUng of condensates 
[13| , explicitly on the particle number N. As examples, 
results for a BEC with A^ = 30 and A^ = 100 bosons are 
presented in Fig. [5][a) and (b), respectively. The decay 
rate in Eq. Q is strongly dominated by the exponential 
term exp(— TV^b/Zi). The action of the bounce trajectory 
tends to zero at the critical value Ocr of the tangent bi- 
furcation [see Fig.|4][b)], and therefore macroscopic quan- 
tum tunneling can play a non-negligible role at scatter- 
ing lengths slightly above that value. The variational 
results in Fig. [5] converge rapidly to those of the grid 
computations with increasing number of coupled Gaus- 
sians. It should be noted that the decay rate obtained 
with a single Gaussian function deviates in the region 
a > a°f — —1.0251 from the converged result by about 
12 orders of magnitude for iV = 30 and by about 20 
orders of magnitude for A^ = 100 particles. 



The decay rates in Fig. [5] are given in particle num- 
ber scaled dimensionless units. For the experiment pro- 
posed by O'Dell et al. |5] experimentally feasible param- 
eters have been found in Ref. [53]. With these values 
the scattering length and the decay rates can be given in 
SI units. For ^''Rb with a.„ = h^/{mu) - 2.5 x 10"^ m 
and tu = 2ma1^/h ^ 27.1 s we present in Table [l] the de- 
cay rates Fg in s"^ at scaled scattering length a = — 1 
corresponding to a^^^^ = — a^Z-ZV^ = —2.78 x 10^''m 
for = 30 and a^^^) = -Uu/N^ = -2.5 x 10"^ m for 
N ~ 100 particles. The values in Table |l] indicate that 
the lifetimes of the self-trapped BEC with reasonable ex- 
perimental parameters can be in the order of seconds at 
scattering lengths close to the tangent bifurcation. The 
strong deviations between the decay rates obtained with 
a single Gaussian function and with improved methods 
are also evident. Thus, the results presented in Fig. [5] 
and Table |l] clearly illustrate that the accurate compu- 
tation of the bounce trajectory is extremely important 
for a reasonable quantitative description of macroscopic 
quantum tunneling. 



V. CONCLUSION AND OUTLOOK 

We have shown that the exact computation of the 
bounce trajectory of a BEC using an extended varia- 
tional approach with coupled Gaussian functions is a full- 
fledged alternative to numerically expensive simulations 
on space-time lattices. Very few (about three to five) 
Gaussians are sufficient to obtain converged results. The 
method here has been developed for a spherically sym- 
metric self-trapped BEC with an attractive l/r interac- 
tion, however, it can straightforwardly also be applied to 
condensates in an external trap or with a different long- 
range interaction. The variational approach may be of 
particular advantage for interactions without spherical 
symmetry, e.g., in dipolar condensates, which require an 
increased dimension of the space-time lattice for numer- 
ical simulations. 

The parameters of the exact bounce trajectory have 
been used as input in the formula ^ derived by Stoof 
|13) . This certainly yields improved, however, still not 
exact values Fq for the decay rate of the BEC by macro- 
scopic quantum tunneling. The reason is that the fluctu- 
ations around the stationary ground state and the time- 
dependent bounce trajectory must be considered, i.e., the 
decay rates Fq in Eq. ([8| must be multiplied with the fluc- 
tuation determinant. The computation of the fluctuation 
determinant and the investigation of macroscopic quan- 
tum tunneling in dipolar condensates is work in progress. 

The discussions in this paper are based on the mean- 
field approach for an effective one-particle condensate 
wave function. For a BEC with attractive interactions 
and not too many particles effects of depletion can, in 
principle, be studied beyond the Gross-Pitaevskii theory. 
In that case the dynamics of the decay should strongly 
depend on the particle number, in particular for low num- 
bers of particles, and the loss of particles will probably 
increase the decay. Such computations and comparisons 
with results obtained within the mean-field approach are 
an interesting task for future work. 
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Appendix A: Gaussian integrals 



For the ansatz with real valued coupled Gaussian func- 



tions given in Eq. ( 21 ) the integrals required in Sec. Ill A 
read 



(gilgk) 
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-7fc! 
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l5/2 
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(Al) 
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FIG. 6. (Color online) Sketch of a potential with an inverted 
harmonic saddle at the origin. 



{9i\'r'^ \9k) = ^TT 
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3 ^7/2 ' 



3/2* 



^kl 



N„ 



-lijkl 



• (A ^3/2 ' 

Ng 



i-J = l ijkl 



-fijki 



i,j=l ^ij^klA^jkl 



(A3) 
(A4) 
(A5) 
(A6) 



A A'^' 
^I'J^kl^ijkl 



with 



Aki = Ak + Ai , 
Aij = Ai + Aj , 



(A7) 



(A8) 
(A9) 



A^Jkl = A, + Aj +Ak + Ai, (AlO) 
7fc/ = 7/c +7; 1 (All) 
Izjki = 7i + 7i + 7fc + 7/ • (A12) 

Appendix B: Calculation of vq 

For the calculation of vq we assume a potential V{x) 
with an inverted harmonic saddle at the origin as illus- 
trated in Fig. |6] The period of a trajectory at energy 
A£; < is written as T{AE) = 2(To + Ti), with Tq the 
time in the region [a:;o,xi] where the harmonic approxi- 
mation is valid (see Fig. IgI . We obtain 



To 



dx 

X 



UJ 







X, 







1 , -1 a;i xo<s^xi 1 2xi 
= — cosh — w — m 

UJ Xq 

T(AE) = 2(To + Ti) « - In f ^e^^' 

UJ \ Xq 



(Bl) 



When the turning point xq = \/—2AEjmuj^ and the 
point xi given by [see Eq. (Ill] 



xi = x(Ti) 



(B2) 



are inserted, Eq. (Bl| yields 



T{AE) 



2^ /2^X 1^/2^1 

UJ \LOXq J LO \—AE 



Wo = 



-AS 
2m ^ 



jT/2 



(B3) 



With T = w = Wo, A£: = - £;mf(/3), and taking 
the limit /? — )■ cx) to approach the bounce trajectory we 
end up with Eq. ( 32 ) . 
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